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The anomalous scaling behavior of the n-th order cor- 
relation functions T n of the Kraichnan model of turbulent 
passive scalar advection is believed to be dominated by the 
homogeneous solutions (zero-modes) of the Kraichnan equa- 
tion B n Tn = 0. Previous analysis found zero-modes in per- 
turbation theory in a small parameter. We present non- 
perturbative analysis of the simplest (non-trivial) case of 
n — 3 and compare the results with the perturbative pre- 
dictions. 



The Kraichnan model of turbulent passive scalar ad- 
vection |jj has attracted enormous attention recently, 
being the first non-trivial model of turbulent statis- 
tics in which the phenomenon of multi-scaling seems un- 
derstandable by analytic methods. The model is for a 
scalar field T(r, t) which satisfies the equation of motion 



dT{r,t) 

at 



+ u(r, t) • VT(r, t) = nV 2 T{r, t) + £(r, t). (1) 



Here £{v,t) is a Gaussian white random force, n is the 
diffusivity and the driving field u(r,t) is chosen to have 
Gaussian statistics, and to be "fastly varying" in the 
sense that its time correlation function is proportional 
to 5(t). The statistical quantities that one is interested 
in are the many point correlation functions 

■F 2 „(ri,r 2 ,...,r 2n ) = ((T(n, t)T(r 2 , t) . . . T{r 2n , <))), (2) 

where double pointed brackets denote an ensemble aver- 
age with respect to a stationary statistics of the forcing 
and the statistics of the velocity field. One of Kraich- 
nan's major results j^j is an exact differential equation 
for this correlation function, 



RHS. 



(3) 



The operator 2? 2ra = 5Z Q >aS Qj 3, and B a are defined by 

B a p ee B(Y a ,Yp) = hij(r a - v^d 2 /dr a ^drp^ , (4) 

where the "eddy-diffusivity" tensor h{j (R) is given by 

hij(R) = h(R)[(( h + d- 1)<% - ChRiRj/R 2 }, 

and h(R) = H(R/C)^ h , < ( h < 2. Here C is some 
characteristic outer scale of the driving velocity field. The 
scaling properties of the scalar depend sensitively on the 



scaling exponent C,h that characterizes the R dependence 
of hij(R) and that can take values in the interval [0,2]. 
Finally, the RHS in Eq. (|J) is known explicitly, but is not 
needed here. The reason is that it was argued that the 
solutions of this equation for n > 1 are dominated by the 
homogeneous solutions ("zero-modes"), in the sense that 
deep in the inertial interval the inhomogeneous solutions 
are negligible compared to the homogeneous one. Also, it 
was claimed that in the inertial interval one can neglect 
the Laplacian operators in Eq.(^), and remain with the 
simpler homogeneous equation B 2 „.F 2n = 0. 

Having exact differential equations for J-% n allowed 
Kraichnan to announce a mechanism for anomalous 
scaling 0]. Assuming that the physical solutions are 
scale invariant one needs to examine the scaling (or ho- 
mogeneity) exponent £ 2 „ of Tin which is defined by 
^2n(Ari, Ar 2 . . . Ar 2n ) = A C2 "J c ' 2 „(ri, r 2 . . . r 2n ) if such 
a solution exists. One expects it to exist in the inertial 
range, i.e. all the separations r,j satisfy 77 <C Tij <C L 
where rj and L are the inner and outer scales respec- 
tively. It is known [Q that for T% such a solution ex- 
ists with £ 2 = 2 — (h- If one solves for these expo- 
nent for n > 1, one can understand, at least in this 
simple model, what are the mechanisms for deviations 
from the predictions of dimensional analysis, with possi- 
ble insight also for the Navier-Stokes problem. In search- 
ing methods for computing these important exponents, 
there emerged two basic strategies. One strategy con- 
sidered the differential equation in the "fully unfused" 
regime in which all the separations between the coordi- 
nates are in the inertial range. In this case even in the 
simplest case of n = 2 the function T\ depends on six 
independent variables (for dimensions d > 2), and one 
faces a formidable analytic difficulty for exact solutions. 
Accordingly, several groups considered perturbative so- 
lutions in some small parameter, like (^] or the inverse 
dimensionality 1/rf 0. The rationale for this approach is 
that at £), = and d — > 00 one expects "simple scaling" 
with £ 2n = n£ 2 . The exponent £4, and later also the set 
C2n , were computed as a function of near these simple 
scaling limits. Another strategy considered the differen- 
tial equation in the "fully fused" regime, in which the 
correlation function degenerates to the structure func- 
tion S 2n (R) = ({[T(r + R) — T(r)] 2n )). In this approach 
there is an enormous simplification in having only one 
variable, but one loses information in the process of fu- 
sion. The lost information was supplemented Q by a yet 
underived conjecture about the properties of conditional 
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averages, leading at the end to a close-form calculation 
of the exponents £ 2n for arbitrary dimension and values 
of Oi ■ The results of the two strategies are not in agree- 
ment. Even though numerical simulations |^| and also 
experiments lend support to the assumption used in 
the second strategy and to its computed values of ( 2n , 
there remains an important mystery as to why the two 
approaches reach such different conclusions. The aim of 
this Letter is to explore non-perturbative calculations of 
the zero modes and their exponents, to shed further light 
on this issue. 

Our strategy is to solve exactly, eigenfunctions in- 
cluded, the homogeneous equation satisfied by the 3'rd 
order correlation function ^(rx, r 2 , r 3 ). Note that in 
Kraichnan's model all the odd-order correlation functions 
J~2n+i are zero because of symmetry under the transfor- 
mation T — > — T. This symmetry disappears for example 
H if the random force £(r,i) is not Gaussian (but 8- 
correlated in time), and in particular if it has a non-zero 
third order correlation 

V 3 ( ri ,r 2 ,r 3 )= /dMi2<£(ri,ti)£(r 2 ,t2)£(r3,0)). (5) 

With such a forcing the third order correlator is non-zero, 
and it satisfies the equation 

B 3 J r 3 {r 1 ,r 2 ,r 3 )=V 3 , B 3 = B 12 + B 13 + B 23 . (6) 

This equation pertains to the inertial interval and ac- 
cordingly we neglected the Laplacian operators. We also 
denoted V 3 = \\m raf} ^o V 3(7*1, r 2 , r 3 ). The solution of 
this equation is a sum of inhomogeneous and homoge- 
neous contributions, and below we study the latter. We 
will focus on scale invariant homogeneous solutions which 
satisfy T${\rx, Ar 2 , Ar 3 ) = A 1 ' 3 J- 3 (ri, r 2 , r 3 ). We refer 
to these as the "zero modes in the scale invariant sec- 
tor" . We note that the scaling exponent of the inhomo- 
geneous scale invariant contribution can be read directly 
from power counting in Eq.(^) (leading to £3 = £2)- Any 
different scaling exponent can arise only from homoge- 
neous solutions that do not need to balance the constant 
RHS. In addition, note that scale- invariant zero- modes 
arise not only due to the omission of the diffusive terms 
from Eq. (^), but also as a result of the omission of the 
boundary conditions for large separation (at the outer 
scale L) . The smooth connection to either small or large 
scales must ruin scale invariance. The scale invariant so- 
lutions of Eq.(||) live in a projective space whose dimen- 
sion is lowered by unity compared to the most general 
form; These solutions do not depend on three separa- 
tions but rather on two dimcnsionlcss variables that are 
identified below. It will be demonstrated how boundary 
conditions arise in this space for which the operator B 3 
is neither positive nor self-adjoint. 

Equation (^) is also invariant under the action of the 
d dimensional rotation group SO(d), and under permu- 
tations of the three coordinates. Here we seek solutions 



in the scalar representation of SO(d), where the solution 
depends on the 3 separations ri 2 , r 23 and r 3 i only. We 
transform coordinates to the variables x\ — |r 2 — r 3 | 2 , 
x 2 = 1 7*3 — J"i| 2 , x 3 = |ri — r 2 | 2 . The triangle inequalities 
in the original space are equivalent to the condition 

2(xix 2 + x 2 x 3 + x 3 xi) >x\+x\+ x\. (7) 

The advantage of the new coordinates is that the inequal- 
ity ^ describes a circular cone in the X\, x 2 , x 3 space 
whose axis is the line x\ = x 2 = x 3 and whose circular 
cross section is tangent to the planes x\ — 0, x 2 — and 
x 3 = 0. This cone can be parameterized by three new 
coordinates s, p, 4>: 

x n = s{l — pcos[0 + (27r/3)n]}, 

< s < 00, < p < 1, < 4> < 2tt. (8) 

The s coordinate measures the overall scale of the triangle 
defined by the original coordinates, and configurations 
of constant p and <j) correspond to similar triangles. The 
p coordinate describes the deviation of the triangle from 
the equilateral configuration (p = 0) up to the physical 
limit of three collinear points attained when p = 1; (f> 
does not have a simple geometric meaning. 

The transformation of the linear operator B 3 to the 
new coordinates is straightforward, and produces a sec- 
ond order linear partial differential operator in the s, p, <p 
variables (the full form of the operator is long and will not 
be given here). The scale invariant solution take on the 
form s^ 3 ' 2 f(p,(j)), and the transformed operator applied 
to this form gives an equation for f(t, 4>) 

4(C3)/(p, 4>) = Hp, + b( P , 4>)dl + c( P , <j>)d p d<, (9) 

+u(p, 0, C3)d p + V(p, (j>, C3)<90 + w(p, (j), C3)]/(p, (t>) = • 

The new operator B 3 depends on £ 3 as a parameter and 
it acts on the unit circle described by the polar p, (f> coor- 
dinates. The circle represents the projective space of the 
physical cone described above. 

The discrete permutation symmetry of the original 
Eq.(||) results in a symmetry of Eq.(||) with respect to 
the 6 element group generated by the transformation 
<f> — > <p + 2tt/3 (cyclic permutation of the coordinates 
in the physical space) and (j) —> —(j) (exchange of coordi- 
nates). This symmetry extends to a full U(l) symmetry 
in the two marginal cases of Qh — and C/i = 2 for which 
all the coefficients in (||) become ^-independent. The 
coefficients in (||) all have a similar structure, and for 
example a(p, </>) reads 

a(p, 4>) = ^[1 - pcos(0 + §7m)] (Ch - 2)/2 5(p, (j) + §7rn) , 

n 

where a(p,<fi) is a low order polynomial in p, cos<^> and 
s'm(j> which vanishes at p = 1, (/> = 0. We see that the 
coefficients are analytic everywhere 
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FIG. 1. The scaling exponent £3 as a functions of Ch found 
as the loci of zeros of the determinant of the matrix B3, for 
d = 2. 

on the circle except at the three points p = 1, <f> = 27rn/3 
where n = 0, 1, 2. These points correspond to the fusion 
of one pair of coordinates, and the coefficients exhibit a 
branch point singularity there. This singularity leads to 
a nontrivial asymptotic behavior of the solutions which 
had been described before in terms of the fusion rules 
Note that for Ch = 2 the singularity disappears 
trivially. For Ch = there is also no singularity since a 
exactly compensates for the inverse power. 

The boundary conditions follow naturally when one 
realizes that B3 is elliptic for points strictly inside the 
physical circle. This is a consequence of the ellipticity of 
the original operator B3. On the other hand B 3 becomes 
singular on the boundary p = 1, where the coefficients 
a(p, 4>) and c{p, <fi) vanish. This singularity reflects the 
fact that this is the boundary of the physical region. It 
follows that B 3 restricted to the boundary becomes a re- 
lation between the function f(p — l,(f>) = g{<j)) an d its 
normal derivative d p f{p,4>)\ P =i = h((f>). The relation is 
bg" +uh + vg' + wg = 0. Solutions of Eq.(||) which do not 
satisfy this boundary condition are singular, with infinite 
p derivatives at p = 1. Such solutions are not physical 
since they involve infinite correlations between the dissi- 
pation (second derivative of the field) and the field itself 
when the geometry becomes collinear, but without fu- 
sion. 

Having a homogeneous equation with homogeneous 
boundary conditions we realize that non-trivial solutions 
are available only when det(i?3) = 0. This determinant 
depends parametrically on (3. Since the operator is de- 
fined on a compact domain we expect the determinant 
to vanish only for discrete values of £3 for any given val- 
ues of Ch and the dimensionality d. We know that there 
always exists a trivial constant solution associated with 
C3 = 0. Our aim is to find the lowest lying positive real 
solutions C3 for which the determinant vanishes. 

FIG. 2. Same as Fig.l, but for d = 3. 

FIG. 3. Same as Fig.l, but for d = 4. 

We approach the problem numerically by discretizing the 
operator B3 including the boundary conditions, and solv- 
ing the analogous problem for the discretized operator. 
Using the symmetry of the problem we restricted the do- 
main to one sixth of the circle, and defined a nine-point 
finite difference scheme for the evaluation of the second 
order derivatives. The discretized boundary conditions 
at p = 1 were achieved with the same scheme. The sym- 
metry implies that the new boundary conditions on the 
lines = 0, 7r/3 are simple Neuman boundary conditions 



d<j,f(p, 4>) = 0. After discretization the problem trans- 
forms to a matrix eigenvalue problem B^ = 0, where 
B3 is a large sparse matrix, whose rank depends on the 
mesh of the discretization, and *f? is the discretized /. 
We used NAG's sparse Gaussian elimination routines to 
find the zeros of det(i? 3 ), and determined the values of 
£3 for these zeros as a function of Ch ■ The results of this 
procedure for space dimensions d — 2,3,4 are presented 
in Figs. 1,2, and 3. 

The various branches shown in Figs. 1-3 can be orga- 
nized on the basis of the perturbation theory of the type 
proposed in |i| near Ch = 0. We performed that type of 
analysis and found that at Ch = the allowed values of 
£3 are organized in two sets, 

C 3 f (m,?i) = 2(3m + 2n), 

C3 (m,n) = -2(d-l + 3m + 2n), (10) 

where n and m are any non- negative integer. The low- 
est lying positive values are 4, 6, 8 etc, whereas for d = 2 
the highest negative value is —2. We see that the non- 
perturbative solution displays in all dimensions a branch 
(dashed line) which begins at Ch = 0, £3 = 4 and ends 
at Ch = 2, £3 = 0. This branch is identical to the lowest 
lying positive branch predicted by the perturbation the- 
ory. We computed the slope of this branch near Ch = in 
perturbation theory, and found that it is 2(2 — d)/(d— 1), 
in agreement with the numerics. Also the slopes of the 
other branches that begin at Ch — were obtained per- 
turbatively and found to agree with the numerics. The 
negative branch (shown only for d — 2) never rises above 
its perturbative limit and is not relevant for the scaling 
behaviour at any value of Ch- Note also that the point 
Ch = 2, £3 = appears to be an accumulation point of 
many branches, and we are not confident that all the 
branches there were identified by our finite discretiza- 
tion scheme. This raises a worry about the availability 
of a smooth perturbative theory around Ch = 2. At least 
we expect such a perturbation theory to be very singu- 
lar. Preliminary analytical work indicates that all the 
branches join the point Ch = 2, £3 = with an infinite 
slope. 

The results of our nonperturbative approach lend sup- 
port to the validity of the perturbative calculations of the 
zero modes of -B4. The disagreement between the scaling 
exponents £4 an d the higher order exponents £„ com- 
puted via the perturbative approach and the predictions 
of the other approach based on the fully fused theory 
cannot be ascribed to a formal failure of the perturba- 
tion theory. There are therefore a few possibilities that 
have to be sorted out by further research: 

(i) The crucial assumption that goes to the fully fused 
approach, which is the linearity of the conditional aver- 
age of the Laplacian of the scalar, is wrong. 

(ii) The computation of the zero modes which is achieved 
by discarding the viscous terms in B n is irrelevant for the 
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physical solution. It is not impossible that the diffusive 
term act as a singular perturbation on some of the scale 
invariant modes. In fact, the operator B n with the vis- 
cous term is positive definite and it has no zero modes. 
That this is a possibility is underlined by recent calcu- 



lations of a shell model of the Kraichnan model |13|, in 
which it was shown that the addition of any minute diffu- 
sivity changes the nature of the zero modes qualitatively, 
(iii) Lastly, and maybe most interestingly, it is possible 
that the physical solution is not scale invariant [jl2|. In 
other words, it is possible that ^(ri, r 2 , r 3 ) is not a ho- 
mogeneous functions with a fixed homogeneity exponent 
£3, but rather (for example), that £3 depends on the ratios 
of the separations (or, in other words, the geometry of 
the triangle defined by the coordinates) . If this were also 
the case for even correlation functions J~2n, this would 
open an exciting route for further research to understand 
how non-scale invariant correlation functions turn, upon 
fusion, to scale invariant structure functions. 

In light of the numerical results of ref. || and the ex- 
perimental results displayed in we tend to doubt 
option (i). If we were to guess at this point we would 
opt for possibility (ii). More work however is needed to 
clarify this important issue beyond doubt. 
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